##################################################################
# Nadiya Kostyuk
# Robustness Checks
# R version: 4.1
# 11/19/25
##################################################################

rm(list=ls())

## Install & load packages (all at once)
list.of.packages <- c('dplyr', 'MASS', 'lmtest', 'stargazer', 'stringr', 
                      'lubridate')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}
lapply(list.of.packages, require, character.only = TRUE)


# setting directory: 
setwd("")


# loading data
load('data/data_final.RData')

# subsetting datasets: 
data_d_contested_80 = data_final %>% 
  filter(territorial_control_80_1dlag == 'Contested')

data_d_no_ruscontrol_80 = data_final %>% 
  filter(territorial_control_80_1dlag != 'Russian occupation')
data_d_no_ruscontrol_70 = data_final %>% 
  filter(territorial_control_70_1dlag != 'Russian occupation')


########################################
# Section 2.3 in the Online Appendix
# Robustness Checks + 
# Alternative Explanations
#########################################
# 
#########################################
# 2.3.1 Alternative territorial control
# 70%
##########################################

mod_list = list()

mod1a = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(war_phases_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_70, 
  family = binomial()
)
summary(mod1a)
mod_list[[1]] = mod1a


mod1aa = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(territorial_control_70_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_70, 
  family = binomial()
)
summary(mod1aa)
mod_list[[2]] = mod1aa


#########################################
# 4.2. Alternative model specifications
# probit + fixed effects
##########################################


mod1b = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(war_phases_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial(link = "probit")
)
summary(mod1b)
mod_list[[3]] = mod1b


mod1bb = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(territorial_control_80_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial(link = "probit")
)
summary(mod1bb)
mod_list[[4]] = mod1bb

# saving model results: 
#saveRDS(mod_list, file = 'analysis/output/tables/robcheck1.RDS')

###########################
# Table 5
# Online Appendix
###########################
stargazer(mod1a, mod1aa, mod1b, mod1bb)


##########################
# fixed-effects: 
mod1c = glm(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(war_phases_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + as.factor(ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial()
)
summary(mod1c)

mod1cc = glm(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(territorial_control_80_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + as.factor(ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial()
)
summary(mod1cc)


############################################
# 4.3. Sensitivty Analysis:
# alt measure of military strategy
############################################

ff = data_final %>%
  mutate(
    rus_war_phases_new = case_when(
      control_ratio_moving_average <= -0.002 ~ 'Rus Maneuver Offensive', 
      control_ratio_moving_average >= 0.002 ~ 'Rus Maneuver Defensive', 
      TRUE ~ war_phases # 'Attrition'
    )
  ) %>% 
  arrange(ADM1_NAME, DATE) %>%
  group_by(ADM1_NAME) %>%
  mutate(rus_war_phases_new_1dlag = lag(rus_war_phases_new, n = 1)) %>%
  ungroup()

# subsetting data:
ff_contested_80 = ff %>% 
  filter(territorial_control_80_1dlag == 'Contested')
ff_80 = ff %>% 
  filter(territorial_control_80_1dlag != 'Russian occupation')


mod_list = list()

mod1d = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(rus_war_phases_new_1dlag) 
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = ff_contested_80, 
  family = binomial()
)
summary(mod1d)
mod_list[[1]] = mod1d

mod1dd = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + territorial_control_80_1dlag
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = ff_80, 
  family = binomial()
)
summary(mod1dd)
mod_list[[2]] = mod1dd



#########################################
# 4.3. Alternative explanations
# air strikes
##########################################
data_d_contested_80 = data_d_contested_80 %>%
  mutate(log_rus_airstrike_t_1dlag = log(rus_airstrike_t_1dlag+1))

data_d_no_ruscontrol_80 = data_d_no_ruscontrol_80 %>%
  mutate(log_rus_airstrike_t_1dlag = log(rus_airstrike_t_1dlag+1))


mod1e = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(war_phases_1dlag)
  + log_rus_airstrike_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_contested_80, 
  family = binomial()
)
summary(mod1e)


mod1ee = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + territorial_control_80_1dlag
  + log_rus_airstrike_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial()
)
summary(mod1ee)
mod_list[[4]] = mod1ee

# saving model results: 
#saveRDS(mod_list, file = 'analysis/output/tables/robcheck2.RDS')

stargazer(mod1d, mod1dd, mod1e, mod1ee)


#########################
# 2.3.5 Section 
# (Online Appendix)
# No lagged DV
##########################

mod_list = list()

mod1f = glmer(
  form = disruption_intentional_3indicators  ~
    as.factor(war_phases_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_contested_80, 
  family = binomial()
)
summary(mod1f)
mod_list[[1]] = mod1f


mod1ff = glmer(
  form = disruption_intentional_3indicators  ~
    as.factor(territorial_control_80_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial()
)
summary(mod1ff)
mod_list[[2]] = mod1ff

# saving model results: 
#saveRDS(mod_list, file = 'analysis/output/tables/robcheck3.RDS')

stargazer(mod1f, mod1ff)


#########################
# 2.3.6 Section 
# (Online Appendix)
# E-value analysis
##########################

install.packages('EValue')
library(EValue)

######################
# Maneuver versus 
# Attrition
######################

mod = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(war_phases_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_contested_80, 
  family = binomial()
)
summary(mod)
# Calculate confidence intervals for the effect estimate
conf_intervals <- exp(confint(mod))

#############################
# Rus Offensive Maneuver
#############################

# Calculate observed association (odds ratio)
#@ one as a time
observed_association <- exp(fixef(mod)["as.factor(war_phases_1dlag)Rus Maneuver Offensive"])
observed_association

# Extract lower and upper bounds of the confidence interval
lower_bound <- conf_intervals["as.factor(war_phases_1dlag)Rus Maneuver Offensive", "2.5 %"]
lower_bound
upper_bound <- conf_intervals["as.factor(war_phases_1dlag)Rus Maneuver Offensive", "97.5 %"]
upper_bound

# Calculate E-value
evalue_results <- EValue::evalues.OR(measure = "OR",
                                    est = observed_association, lo = lower_bound, hi = upper_bound, rare = 1)
evalue_results


#############################
# Rus Defensive Maneuver
#############################

# Calculate observed association (odds ratio)
#@ one as a time
observed_association_def <- exp(fixef(mod)["as.factor(war_phases_1dlag)Rus Maneuver Defensive"])
observed_association_def

# Extract lower and upper bounds of the confidence interval
lower_bound_def <- conf_intervals["as.factor(war_phases_1dlag)Rus Maneuver Defensive", "2.5 %"]
lower_bound_def
upper_bound_def <- conf_intervals["as.factor(war_phases_1dlag)Rus Maneuver Defensive", "97.5 %"]
upper_bound_def

# Calculate E-value
evalue_results_def <- EValue::evalues.OR(measure = "OR",
                                     est = observed_association_def, 
                                     lo = lower_bound_def, hi = upper_bound_def, rare = 1)
evalue_results_def


#########################################
# Punishment Hypothesis
#########################################
mod = glmer(
  form = disruption_intentional_3indicators  ~
    disruption_intentional_3indicators_1dlag
  + as.factor(territorial_control_80_1dlag)
  + log_rus_mil_t_1dlag
  + holiday_dummy
  + as.factor(WEEKEND)
  + (1|ADM1_NAME), 
  data = data_d_no_ruscontrol_80, 
  family = binomial()
)
summary(mod)
# Calculate confidence intervals for the effect estimate
conf_intervals <- exp(confint(mod))

# Calculate observed association (odds ratio)
#@ one as a time
observed_association_tc <- exp(fixef(mod)["as.factor(territorial_control_80_1dlag)Ukrainian control"])
observed_association_tc

# Extract lower and upper bounds of the confidence interval
lower_bound_tc <- conf_intervals["as.factor(territorial_control_80_1dlag)Ukrainian control", "2.5 %"]
lower_bound_tc
upper_bound_tc <- conf_intervals["as.factor(territorial_control_80_1dlag)Ukrainian control", "97.5 %"]
upper_bound_tc

# Calculate E-value
evalue_results_tc <- EValue::evalues.OR(measure = "OR",
                                        est = observed_association_tc, 
                                        lo = lower_bound_tc, hi = upper_bound_tc, 
                                        rare = 1)
evalue_results_tc



###############################
# 2.4. Reporting Biases
###############################

###########################
# 2.4.2. Visualization
###########################


data_final = data_final %>%
  mutate(
    outage_cause = ifelse(disruption_intentional_3indicators == 0, 'no outage', outage_cause), 
    cyber_outage_cause = ifelse(outage_cause == 'cyber'|outage_cause == 're-route', 1, 0), 
    military_outage_cause = ifelse(outage_cause == 'military'|outage_cause == 'equipment dismantling', 1, 0)
    ) 
data_summary = data_final %>%
  filter(disruption_intentional_3indicators != 0 | cyber_outage_cause != 0 | military_outage_cause != 0) %>%
  group_by(war_phase_names, cyber_outage_cause, military_outage_cause) %>%
  summarize(total_outages = n())
data_summary

data_long <- data_summary %>%
  pivot_longer(
    cols = c(cyber_outage_cause, military_outage_cause),
    names_to = "outage_type",
    values_to = "outage_flag"
  ) %>%
  filter(outage_flag == 1) %>%  
  dplyr::select(war_phase_names, outage_type, total_outages)
print(data_long)


# Plotting
data_long %>%
  # updating phases with date: 
  mutate(
    war_phase_names = case_when(
      war_phase_names == 'Ph1: Rus Maneuver Offensive: Initial Invasion' ~ 'Ph1: Offensive Maneuver (R) \n (02/24/2022-04/07/2022)',
      war_phase_names == 'Ph2: Attrition: Southeastern front' ~ 'Ph2: Attrition \n (04/08/2022-08/28/2022)', 
      war_phase_names == 'Ph3: Rus Maneuver Defensive: Ukrainian counteroffensive' ~ 'Ph3: Defensive Maneuver (R) \n (08/29/2022-11/11/2022)', 
      war_phase_names == 'Ph4: Attrition: Second stalemate' ~ 'Ph4: Attrition \n (11/12/2022-06/07/2023)', 
      war_phase_names == 'Ph5: Rus Maneuver Defensive: Ukrainian counteroffensive' ~ 'Ph5: Defensive Maneuver (R) \n (06/08/2023-11/30/2023)', 
      war_phase_names == 'Ph6: Attrition: Winter campaigns' ~ 'Ph6: Attrition \n (12/01/2023-12/31/2023)'
    )
  ) %>%
ggplot(aes(x = war_phase_names, y = total_outages, fill = outage_type)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    x = "", 
    y = "Number of Outages per Region-Day", 
    #title = "Total Number of Outages in the Ukraine Conflict by Cause and War Phase",
    fill = "Cyber Outage"
  ) +
  scale_fill_manual(values = c("gray", "black"), labels = c("Cyber Outage (R)", "Military Outage (R)")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), legend.title = element_blank())  # Rotate phase names 45 degrees
#ggsave(file = 'analysis/output/figures/outage_cause.pdf', width = 10, height = 6)




